function [f] = ObjValue(x)
%ObjValue: Objective Function
%% ALL GLOBAL OPTIMA
%% one-dimensional equal minima
% f = 1-sin(5.*pi.*x).^6;
% % d = 1;
% % x_bound = [0,1];
% % OptSet.sol = [0.1;0.3;0.5;0.7;0.9];
% % OptSet.fval = 0;
% % OptSet.radius = 0.015;
% % precision_init = 0.04; % 0.2
% % precision = [0.023,0.0075,0.0023,0.00075];
% % MaximalBudget = 50000;
%% Himmelblau's function
% f = (x(:,1).^2+x(:,2)-11).^2+(x(:,1)+x(:,2).^2-7).^2;
% % d = 2;
% % x_bound = repmat([-6,6],[d,1]);
% % OptSet.sol = [3,2;
% %     -2.805118,3.131312;
% %     -3.779310,-3.283186;
% %     3.584428,-1.848126];
% % OptSet.fval = 0;
% % OptSet.radius = 0.1;
% % precision_init = 0.4; % 3.848
% % precision = [0.074, 0.024, 0.0074, 0.0024];
% % MaximalBudget = 50000;
%% Six-hump camel back
% f = 4.*((4-2.1.*(x(:,1).^2)+x(:,1).^4./3).*(x(:,1).^2)+x(:,1).*x(:,2)+(4.*(x(:,2).^2)-4).*(x(:,2).^2));
% % d = 2;
% % x_bound = [-1.9,1.9;-1.1,1.1];
% % OptSet.fval = -4.12651;
% % OptSet.sol = [0.0898207939233943,-0.712627558313027;
% %     -0.0898442050683038,0.712633095015521];
% % OptSet.radius = 0.16;
% % precision_init = 0.15; % 1.4
% % precision = [0.085, 0.028, 0.0085, 0.0028];
% % MaximalBudget = 50000;
%% Modified Rastrigin function
% k = [3,4];
% xn = size(x,1);
% f = sum(10+9.*cos(2.*pi.*repmat(k,[xn,1]).*x),2);
% % d = d;
% % k = [2,2,2,2,2];
% % x_bound = repmat([0,1],[d,1]);
% % solutionId = fullfact(k);
% % OptSet.sol = (2.*solutionId-1)./2./repmat(k,[prod(k),1]);
% % OptSet.fval = d;
% % k=[3,4]: OptSet.radius = 0.01;
% % precision_init = 0.04; % 0.25
% % precision = [0.0096, 0.0030, 0.00096, 0.0003];
% % MaximalBudget = 200000;
% % k=[3,3,3]:  OptSet.radius = 0.01;
% % precision_init = 0.04; % 0.33
% % precision = [0.009, 0.003, 0.0009, 0.0003];
% % MaximalBudget = 400000;
% % k=[2,2,2,2,2]:  OptSet.radius = 0.02;
% % precision_init = 0.07; % 0.5
% % precision = [0.011, 0.0034, 0.0011, 0.00034];
% % MaximalBudget = 400000;
%% Vincent function
% f = 1-mean(sin(10.*log(x)),2);
% % d = d;
% % x_bound = repmat([0.25,10],[d,1]);
% % solutionId = fullfact(ones(1,d).*6);
% % OptSet.sol = exp((4.*solutionId-11)./20.*pi);
% % OptSet.fval = 0;
% % OptSet.radius = 0.03;
% % precision_init = 0.08; % 0.2912
% % precision = [0.03, 0.0096, 0.003, 0.00096];
% % 2D: MaximalBudget = 200000;
% % 3D: MaximalBudget = 400000;
%% Shubert function
% [xn,d] = size(x);
% f = zeros(xn,1);
% a = repmat((1:5)',[1,d]);
% for i = 1:xn 
%     f(i) = prod(sum(a.*cos(repmat(x(i,:),[5,1]).*(a+1)+a),1));
% end
% % d = d;
% % x_bound = repmat([-10,10],[d,1]);
% % OptSet.radius = 0.01;
% % precision_init = 0.16; % 0.6248
% % 2D: OptSet.fval = -186.7309;
% % OptSet.sol = [-7.70831,-7.08351; -7.70831,-0.80032; -7.70831,5.48286; -7.08351,-7.70831;
% %     -7.08351,-1.42513; -7.08351,4.85806; -1.42513,-7.08351; -1.42513,-0.80032;
% %     -1.42513,5.48286; -0.80032,-7.70831; -0.80032,-1.42513; -0.80032,4.85806;
% %     4.85806,-7.08351; 4.85806,-0.80032; 4.85806,5.48286; 5.48286,-7.70831;
% %     5.48286,-1.42513; 5.48286,4.85806];
% % precision = [0.0096, 0.003, 0.00096, 0.0003];
% % MaximalBudget = 200000;
% % 3D: OptSet.fval = -2709.0935;
% % OptSet.sol = [-7.70831,-7.08351,-7.08351; -7.70831,-7.08351,-0.80032; -7.70831,-7.08351,5.48286;
% % -7.70831,-0.80032,-7.08351; -7.70831,-0.80032,-0.80032; -7.70831,-0.80032,5.48286;
% % -7.70831,5.48286,-7.08351; -7.70831,5.48286,-0.80032; -7.70831,5.48286,5.48286;
% % -7.08351,-7.70831,-7.08351; -7.08351,-7.70831,-0.80032; -7.08351,-7.70831,5.48286;
% % -7.08351,-7.08351,-7.70831; -7.08351,-7.08351,-1.42513; -7.08351,-7.08351,4.85806;
% % -7.08351,-1.42513,-7.08351; -7.08351,-1.42513,-0.80032; -7.08351,-1.42513,5.48286;
% % -7.08351,-0.80032,-7.70831; -7.08351,-0.80032,-1.42513; -7.08351,-0.80032,4.85806;
% % -7.08351,4.85806,-7.08351; -7.08351,4.85806,-0.80032; -7.08351,4.85806,5.48286;
% % -7.08351,5.48286,-7.70831; -7.08351,5.48286,-1.42513; -7.08351,5.48286,4.85806;
% % -1.42513,-7.08351,-7.08351; -1.42513,-7.08351,-0.80032; -1.42513,-7.08351,5.48286;
% % -1.42513,-0.80032,-7.08351; -1.42513,-0.80032,-0.80032; -1.42513,-0.80032,5.48286;
% % -1.42513,5.48286,-7.08351; -1.42513,5.48286,-0.80032; -1.42513,5.48286,5.48286;
% % -0.80032,-7.70831,-7.08351; -0.80032,-7.70831,-0.80032; -0.80032,-7.70831,5.48286;
% % -0.80032,-7.08351,-7.70831; -0.80032,-7.08351,-1.42513; -0.80032,-7.08351,4.85806;
% % -0.80032,-1.42513,-7.08351; -0.80032,-1.42513,-0.80032; -0.80032,-1.42513,5.48286;
% % -0.80032,-0.80032,-7.70831; -0.80032,-0.80032,-1.42513; -0.80032,-0.80032,4.85806;
% % -0.80032,4.85806,-7.08351; -0.80032,4.85806,-0.80032; -0.80032,4.85806,5.48286;
% % -0.80032,5.48286,-7.70831; -0.80032,5.48286,-1.42513; -0.80032,5.48286,4.85806;
% % 4.85806,-7.08351,-7.08351; 4.85806,-7.08351,-0.80032; 4.85806,-7.08351,5.48286;
% % 4.85806,-0.80032,-7.08351; 4.85806,-0.80032,-0.80032; 4.85806,-0.80032,5.48286;
% % 4.85806,5.48286,-7.08351; 4.85806,5.48286,-0.80032; 4.85806,5.48286,5.48286;
% % 5.48286,-7.70831,-7.08351; 5.48286,-7.70831,-0.80032; 5.48286,-7.70831,5.48286;
% % 5.48286,-7.08351,-7.70831; 5.48286,-7.08351,-1.42513; 5.48286,-7.08351,4.85806;
% % 5.48286,-1.42513,-7.08351; 5.48286,-1.42513,-0.80032; 5.48286,-1.42513,5.48286;
% % 5.48286,-0.80032,-7.70831; 5.48286,-0.80032,-1.42513; 5.48286,-0.80032,4.85806;
% % 5.48286,4.85806,-7.08351; 5.48286,4.85806,-0.80032; 5.48286,4.85806,5.48286;
% % 5.48286,5.48286,-7.70831; 5.48286,5.48286,-1.42513; 5.48286,5.48286,4.85806];
% % precision = [0.002, 0.00064, 0.0002, 0.000064];
% % MaximalBudget = 400000;
%% Conposition function
[xn,d] = size(x);
global o M
% 1: Sphere; 2: Grienwank; 3: Rastrigin; 4: Weierstrass; 5:EF8F2
%% Composition functions 1
% fun_type = [2,2,4,4,1,1];
% sigma = [1,1,1,1,1,1];
% lambda = [1,1,8,8,1/5,1/5];
% % d = 2;
% % x_bound = repmat([-5,5],[d,1]);
% % global o M
% % load('composition_data.mat','o')
% % o = o(1:6,1:d);
% % M = repmat(diag(ones(1,d)),[1,1,6]);
% % OptSet.radius = 0.01;
% % OptSet.fval = 0;
% % OptSet.sol = o;
% % precision_init = 0.16; % 2.76
% % precision = [5E-7,1.6E-8,1.1E-9,3.2E-10];
% % MaximalBudget = 200000;
%% Composition functions 2
% fun_type = [3,3,4,4,2,2,1,1];
% sigma = [1,1,1,1,1,1,1,1];
% lambda = [1,1,10,10,1/10,1/10,1/7,1/7];
% % d = 2;
% % x_bound = repmat([-5,5],[d,1]);
% % global o M
% % load('composition_data.mat','o')
% % o = o(1:8,1:d);
% % M = repmat(diag(ones(1,d)),[1,1,8]);
% % OptSet.radius = 0.01;
% % OptSet.fval = 0;
% % OptSet.sol = o;
% % precision_init = 0.16; % 2.76
% % precision = [1.6E-6,5.6E-8,1.8E-9,5.5E-10];
% % MaximalBudget = 200000;
%% Composition functions 3
fun_type = [5,5,4,4,2,2];
sigma = [1,1,2,2,2,2];
lambda = [1/4,1/10,2,1,2,5];
% d = 2;
% x_bound = repmat([-5,5],[d,1]);
% global o M
% load('composition_data.mat','o','M_CF3_2D')
% o = o(1:6,1:d);
% M = M_CF3_2D;
% OptSet.radius = 0.01;
% OptSet.fval = 0;
% OptSet.sol = o;
% 2D: precision_init = 0.16; % 2.76
% precision = [9.5E-8,2E-9,1.5E-10,4.3E-11];
% MaximalBudget = 200000;
% 3D: precision_init = 0.16; % 2.76
% precision = [9E-8,2.5E-9,1.4E-10,3.8E-11];
% MaximalBudget = 400000;
% 5D: precision_init = 0.16; % 2.76
% precision = [1.1E-7,3E-9,1.5E-10,4.2E-11];
% MaximalBudget = 400000;
% 10D: precision_init = 0.16; % 2.76
% precision = [3.8E-8,1.2E-9,1.1E-10,3.2E-11];
% MaximalBudget = 400000;
%% Composition functions 4
% fun_type = [3,3,5,5,4,4,2,2];
% sigma = [1,1,1,1,1,2,2,2];
% lambda = [4,1,4,1,1/10,1/5,1/10,1/40];
% % d = 3;
% % x_bound = repmat([-5,5],[d,1]);
% % global o M
% % load('composition_data.mat','o','M_CF4_3D')
% % o = o(1:8,1:d);
% % M = M_CF4_3D;
% % OptSet.radius = 0.01;
% % OptSet.fval = 0;
% % OptSet.sol = o;
% % MaximalBudget = 400000;
% % 3D: precision_init = 0.2; % 2.76
% % precision = [0,0,0,0];
%%
n = length(fun_type);
bias = zeros(1,n);
% weights
w = zeros(xn,n);
for i = 1:xn
    w(i,:) = (exp(-sum((repmat(x(i,:),[n,1])-o).^2,2)./2./d./(sigma.^2)'))';
end
[maxw,wid] = max(w,[],2);
w = w.*(1-repmat(maxw,[1,n]).^10);
for i = 1:xn
    w(i,wid(i)) = maxw(i);
end
w = w./repmat(sum(w,2),[1,n]);
% objective function
fbasic = zeros(xn,n);
for i = 1:n
    fbasic0 = fcom(fun_type(i),(x-repmat(o(i,:),[xn,1]))./lambda(i)*M(:,:,i));
    f_max = abs(fcom(fun_type(i),(ones(1,d).*5)./lambda(i)*M(:,:,i)));
    fbasic0 = fbasic0.*2000./f_max;
    fbasic(:,i) = (fbasic0 + bias(i));
end
f = sum(w.*fbasic,2);

%% MANY LOCAL OPTIMA
%% one-dimensional increasing minima
% f = 1-exp(-2.*log(2).*((x-0.08)./0.854).^2).*(sin(5.*pi.*(x.^(3/4)-0.05)).^6);
% % d = 1;
% % x_bound = repmat([0.02,1],[d,1]);
% % OptSet.fval = [1.84754329279890e-07; 0.0513106897498393;
% %     0.229184764278421; 0.495888491374641; 0.748389918860396];
% % OptSet.sol = [0.0796967889120553; 0.246276975868616;
% %     0.449497982732034; 0.679163842395897; 0.930151599444996];
% % OptSet.radius = 0.01;
% % precision_init = 0.02; % 0.2
% % precision = [0.017,0.0052,0.0017,0.00052];
% % MaximalBudget = 1000;
%% Six-hump camel back
% f = 4.*((4-2.1.*(x(:,1).^2)+x(:,1).^4./3).*(x(:,1).^2)+x(:,1).*x(:,2)+(4.*(x(:,2).^2)-4).*(x(:,2).^2));
% % d = 2;
% % x_bound = [-1.9,1.9;-1.1,1.1];
% % OptSet.fval = [-4.12651; -4.12651; -0.86185; -0.86185];
% % OptSet.sol = [0.0898207939233943,-0.712627558313027;
% %     -0.0898442050683038,0.712633095015521;
% %     1.70356679399635,-0.796073100464659;
% %     -1.70353035333493,0.796131617535420];
% % OptSet.radius = 0.1;
% % precision_init = 0.1; % 1.4
% % precision = [0.066, 0.021, 0.0066, 0.002];
% % MaximalBudget = 30000;
%% Rastrigin function
% d = size(x,2);
% A = 10;
% f = A*d+sum(x.^2-A.*cos(2.*pi.*x),2);
% % d = d;
% % x_bound = repmat([-5.12,5.12],[d,1]);
% % LevelNum = 4;
% % level = [0, 0.9949586377, 1.98991223375, 2.984855700];
% % solutionId = fullfact(ones(1,d)*(LevelNum*2-1))-LevelNum;
% % OptSet.sol = level(abs(solutionId)+1);
% % OptSet.sol(solutionId<0) = -OptSet.sol(solutionId<0);
% % OptSet.fval = ObjValue(OptSet.sol);
% % [OptSet.fval,rank] = sort(OptSet.fval);
% % OptSet.sol = OptSet.sol(rank,:);
% % OptSet.radius = 0.1;
% % precision_init = 0.1; % 1
% % 2D: precision = [0.03, 0.01, 0.003, 0.001];
% % MaximalBudget = 30000;
% % 3D: precision = [0.026, 0.008, 0.0026, 0.0008];
% % MaximalBudget = 100000;

%% CONTINUOUS OPTIMA PROBLEMS
%% Schaffer's function
% f = 0.5+(sin(sqrt(sum(x.^2,2))).^2-0.5)./((1+0.001.*(sum(x.^2,2))).^2);
% % d = 2;
% % x_bound = repmat([-10,10],[d,1]);
% % precision_init = 0.4; % 1
% % precision = [0.4, 0.4, 0.2, 0.08];
% % MaximalBudget = 100000;
%% Michalewicz
% f = -sin(x(:,1)).*(sin(x(:,1).^2./pi).^20)+sin(x(:,2)).*(sin(x(:,2).^2.*2./pi).^20);
% % d = 2;
% % x_bound = repmat([0,pi],[d,1]);

%% OTHERS
%% Shekel's foxholes
% a = 16.*(mod(0:24,5)-2);
% b = 16.*(floor((0:24)./5)-2);
% xn = size(x,1);
% f = zeros(xn,1);
% for i = 1:xn
%     temp = 1+(0:24)+(x(i,1)-a).^6+(x(i,2)-b).^6;
%     f(i) = 1./(0.002+sum(1./temp));
% end
% % d = 2;
% % x_bound = repmat([-65.536,65.535],[d,1])
%% Griewank function
% d = size(x,2);
% f = 1+1/4000.*sum(x.^2,2)-prod(cos(x./repmat(sqrt(1:d),[size(x,1),1])),2);
% % d = 2;
% % x_bound = repmat([-5,5],[d,1]);
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% basic functions for composition function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function fcom = fcom(id,x)
%   id: 1: Sphere; 2: Grienwank; 3: Rastrigin; 4: Weierstrass; 5:EF8F2
        xn_f = size(x,1);
        if id == 1
            fcom = sum(x.^2,2);
        elseif id == 2
            fcom = sum(x.^2,2)./4000-prod(cos(x./repmat(sqrt(1:d),[xn_f,1])),2)+1;
        elseif id == 3
            fcom = sum(x.^2-10.*cos(2*pi.*x)+10,2);
        elseif id==4
            fcom = zeros(xn_f,1);
            kk = (0:20)';
            aa = 0.5.^kk;
            bb = 3.^kk;
            fw0 = d*sum(aa.*cos(2*pi.*bb.*0.5));
            aa = repmat(aa,[1,d]);
            bb = repmat(bb,[1,d]);
            for ii = 1:xn_f
                fcom(ii) = sum(sum(aa.*cos(2*pi.*bb.*(repmat(x(ii,:),[21,1])+0.5))))-fw0;
            end
        else
            x = x+1; % (1,...,1) is the minima
            F2x = (x.^2-x(:,[2:d,1])).^2.*100+(x-1).^2;
            fcom = sum(F2x.^2./4000-cos(F2x)+1,2);
        end
%         fcom = -fcom;
    end

end

